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MULTILEVEL AND MULTISCALE SCHEMES FOR FRACTIONAL PARTIAL 

DIFFERENTIAL EQUATIONS * * 

Zhijiang Zhang^ and Weihua Deng^ 


Abstract. The wavelet numerical methods for the classic PDEs have been well developed, but they 
are still not discussed for the fractional PDEs. This paper focuses on investigating the applications of 
wavelet bases to numerically solving fractional PDEs and digging out the potential benefits of wavelet 
methods comparing with other numerical methods, especially in the aspects of realizing precondition¬ 
ing, adaptivity, and keeping the Toeplitz structure. More specifically, the contributions of this paper 
are as follows: 1. the techniques of efficiently generating stiffness matrix with computational cost 
0 ( 2 '') are provided for first, second, and any order bases; 2. theoretically and numerically discuss the 
effective multilevel preconditioner for time-independent equation and multiresolution multigrid method 
for time-dependent equation, respectively; 3. the wavelet multiscale adaptivity is experimentally dis¬ 
cussed and numerically applied to solve the time-dependent (independent) equations. In fact, having 
reliable, simple, and local regularity indicators is the striking benefit of the wavelet in adaptively solving 
fractional PDEs. 

1991 Mathematics Subject Classification. 35R11, 65T60, 65F08, 65M55. 


1. Introdugtion 

The continuous time random walk (CTRW), a fundamental model in statistic physics, is a stochastic process 
with arbitrary distributions of jump lengths and waiting times. When the jump length and/or waiting time 
distribution(s) are/is power law and the second order moment of jump lengths and/or the first order moment 
of waiting times are/is divergent, the CTRW describes the anomalous diffusion, i.e., the super and sub diffusive 
cases; and its Fokker-Planck equation has space and/or time fractional derivative(s) [22]. It can be noted that 
the corresponding fractional PDEs are essentially dealing with the multiscale phenomena; and generally the 
fractional PDEs have weaker regularity at the area close to boundary and initial time. Besides anomalous 
diffusion, the fractional models are also used to characterize the memory and hereditary properties inherent in 
various materials and processes and, recently, much more scientific applications are found in a variety of fields; 
see, e.g., [T^l^[271139] and the references therein. 

The obtained analytical solutions of fractional PDEs are usually in the form of transcendental functions or 
infinite series; and in much more cases, the analytical solutions are not available. Then the approximation 
and numerical techniques for solving the fractional PDEs become essential and have been developed very fast 
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recently, such as, the finite difference method [sinniisniiMiiii], the finite element method [SHniisliss], and 
the spectral method [siiniiiQ]. But the computational expenses and nonuniform regularity are still the big 
challenges that one faces in numerically solving the fractional PDEs, owing to the nonlocality and potential 
multiscale characteristics of the fractional derivatives; and basing on the preconditioning, adaptivity, and fast 
transform techniques to develop high efficient methods seems to be a new trend. The preconditioning techniques 
are discussed in P1I57]. where the Krylov subspace projection is their common theme. Fast transform method 
and multigrid method are provided in |34) and |24) . respectively. 

So far, there seems to be very limited works [T3l[28l[29l[33] to solve the fractional PDEs or ODEs by wavelet, 
although the wavelet numerical methods for classical PDEs or ODEs have been well developed [3131] . The goal 
of this paper is to dig out the potential advantages of wavelets in treating the fractional operators, including 
preconditioning, multigrid, adaptivity, and keeping the quasi-Toeplitz structure for arbitrary order wavelet 
bases. More concretely, the clearly obtained benefits of wavelets for fractional operator consist of the following: 
1) stiffness matrix of fractional operator is Toeplitz for scaling bases because of their shift-invariant property 
(it is not always true for the familiar finite element bases, such as the quadratic or cubic element) and a simple 
diagonal scaling usually produces a good preconditioner; 2 ) multiscale coefficients indicate the local regularity, 
and they can be used as the indicator (local posteriori error estimate seems hard to be obtained for the adaptive 
finite element method because of the global property of the operator) in the adaptive mesh refinement for 
controlling the entire computational process and increasing the efficiency, i.e., one only needs to make the local 
refinement on the subdomain where the wavelet coefficients are larger compared with those of other places. For 
avoiding all non indispensable complications, we present the main ideas and techniques in their simplest form 
and restrict ourselves to the following homogeneously space fractional PDE [lOl 1^132] : 

qut + Am = / on D, (1.1) 

where D = (0,1), 9 = 0 or 1; and A is a (2 — /?)-th (0 < /3 < 1) order differential operator 

Am := -KjjD (p oD~^ -|- (1 - p) xDi^^ Du (1.2) 

with Kjs > 0 being the generalized diffusivity; 0 < p < 1 , D represents a single spatial derivative; and 

xDi^ are the left and right fractional integral operators [ 22 ], being, respectively, defined as 

oD-^m: = {x- sf-^u{s)ds, (1.3) 

xD~^u: = J is-xf-^u{s)ds. (1.4) 

When g = 1, one gets the fractional initial boundary value problem (IBVP) with an additional initial condition 
u{x, 0) = g{x); and when p = 0, it is the fractional boundary value problem (BVP), which can also be regarded 
as the steady state equation of the associated IBVP. Considering the homogeneous boundary condition and 
using integration by parts, one can easily get DqDy^Du = qD‘^~^u and D^D^^Du = xD1~^u. Then one can 
reduce the model (HU to a more familiar form, and a basic theoretical framework for its variational solution 
has been presented in m-, this will enable us to put focus on the wavelet numerical methods themselves. 

This paper is organized as follows. In Section 2, we give a brief recall to the spline scaling and wavelet 
functions. They have the closed-form expression, which is of course attractive for the fractional operators. In 
Section 3, we study the computational formulation with respect to the uniform grids. We first discuss the 
effective way to construct the algebraic system, and then derive the multilevel wavelet preconditioning and the 
multiresolution multigrid schemes (MMG) for solving the BVP and IBVP, respectively. In Section 4, we give 
some heuristically adaptive algorithms and show how singularities can be easily detected by wavelet, and put 
our attention on its efficiency by proposing and testing the adaptive algorithms that concentrate the degrees 
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of freedom in the neighborhood of near singularities. The numerical results are shown in Section 5 and we 
conclude the paper with some remarks in the last section. 


2. Preliminaries 

In this section, we collect/present some essential properties of the scaling functions and wavelets to make the 
paper self-contained and more readable. For the details, refer to HHS] and mu- First, we give the definitions 
of the fractional Sobolev spaces used in this paper. For any s > 0, let be the Sobolev space of order s 

on R, and the space of the restriction of the functions from H®(R). More specihcally. 



nffR) = |u(x) G L‘^{R) |u|7(r) < ooj 

(2.1) 

endowed with the 

seminorm 



l^lw'»(R) = / \-^[u\{uj)f div 

(2.2) 

and the norm 


ll«ll«RR)= / {l + \LOnmu]iuj)fduJ^ / {l+\uj\^)'^mu]{Lv)\^dcV, 

Jr Jr 

(2.3) 


'H®(P) = {u € 311 G 7d®(R) such that ft o = u} 

(2.4) 

endowed with 


|M|w'>(n) = |m|w«(r) and |lM||7(n) = l|w|lz, 2 (n) + |w|wRn)) 

(2.5) 


where ^[u] denotes the Fourier transform of u. And Ho(^) defined as the closure of C^(fi) w.r.t. || • ||?i»(a). 

Let [xo,..., Xd]f denote the d-th order divided difference of / at the points xq, ..., Xd, := (max{0, t})^ d > 
2; and choose the Schoenberg sequence of knots 


P := {0,...,0,2-J,2 X 2-J,3 X 2-^, 

...,l-2-^^ 


(2.6) 

d 


d 


to define the scaling function sets /c G = {l,... 

.,2^ -kd-3} } 

with 


(t>j,k{x) := 2= (t7d+i ~ ^fc+i)[^fc+i> ■ 

• ■ ’ ^i+d+l](^ ~ 


(2.7) 


which is the scaled B-Splines [55]. Then the sequence Sj = span{$j} forms a multiresolution analysis (MRA) 
of L 2 {I), where / = (0,1) . The system is uniformly local and locally finite, i.e., diam(supp(/j,fc) ~ and 
ff{4'j,k ■ supp</j_fc n supp(/)jy} ^ 1; it forms a stable Riesz basis of Sj, i.e.. 


c$||cj||z2(a^) < ^ < C'$||cj||/2(Aj); 

" fcsA,- 


and Sj satisfies the Jackson and Bernstein estimates, i.e.. 


Jnf ||x - ^2 lkll«^(/) Vx e n^{I) 




( 2 . 8 ) 


< 


\/vj € Sj^ 0 < 5 < 7 , 


(2.9) 

( 2 . 10 ) 
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where 7 := sup{z/ G IR : Vj G Mvj G Sj} and by ^ ^ i? we mean that A can be bounded by a multiple of 

i?, independent of the parameters they may depend on. 

Since is a Riesz basis of Sj^ there exists a dual MRA sequence Sj = span{<I)j}, which also forms a MRA 
of L 2 {fi). And one can define the biorthogonal projector: 

P, : L 2 {n) ^ S„ P,v := ^ (2.11) 

keAj 


where (•, •) denotes the inner product. Then — Pj , ciiici for 0 ^ s ^ 1 j 




0 < s < 7 < d. 


( 2 . 12 ) 


One can also construct the interval biorthogonal wavelet sets 4'^ = {ipj^k^k G V^} and 4'j = G V^}; it 

holds the norm equivalence, i.e., there exist a,a > 0 such that for the Sobolev space 7d®(0), 


3>Jo—i feeVj 


2 


Jo —1 fcGVj 


V s G {-a, a), 


(2.13) 


where ■= ^Jo,fc,Vjo_i := Ajg,djg-i^k ■= cjg^k', Jo denotes the lowest level. It also means that 

U^jo-i is a Riesz basis of 'Ho(^)- 

Moreover, denote Wj = spanlTj}. Then the operator Qj := Pj+i — Pj is a projection onto the space Wj, 
having the representation 

QjV= 'fl^3,k- (2-14) 

k^V j 

And for 0 < 7 < d, there exists 


(/> ^3,k) 


' inf \\f-p\\ 

P&Pd-i 


L2(supp'03,fc) 




’3,k 


A 2-f'^ 


f(7) 


L2(suppi/;3,fc) 


(2.15) 


This shows that the wavelet coefficients are small provided that the function is locally smooth, which is the 
foundation to design wavelet adaptive algorithms. 

Note that based on the scaling function sets some other special wavelets can also be constructed. If 
one demands Wj = Wj and Sj = Sj, the semiorthogonal wavelets can be obtained [7]; if one chooses 
d = 2, tpix) = -^cj)i^i{x),"i’j = {V’j.fe = 2i'il){2^x — k),k G Vj}, then the interpolation wavelet 
is obtained, and it also satisfies the norm equivalence for s G (1, |) (Page 605 of [5]); and if d = 4, let 
= {<(>j,fc, k G Aj/{1, 2f + 1}}, and define 4’{x), (jibix), V'(a^) and ipbix) by 


4>h{.x) = -x+-—xl +-{x-l)l--{x-2)l + ^— - - 

-ipix) = P(l){2x) + (l){2x-l)-^(j){2x-2), 

iib{x) = (j)b{2x) - ^4>(.‘^x). 


(2.16) 

(2.17) 

(2.18) 
(2.19) 


One has the semi-interpolation spline wavelet 'i'j = {tpb{‘2^x),tpj^k |o<fe<( 23 - 3 )i '0b(2f(l — J^))}, which satishes 
the so-called point value vanishing property, and the wavelet expansion coefficients also indicate the regularity 
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of the approximation function [5]. In practice, the base functions ipj.ii') are usually added for 

removing the limitation that the first order derivative of the (to be approximated) function at the boundary 
needs to be zero. 

Finally, we point out that there exists the refinement relations 

ihj = 0, = (2.20) 

And the space Sj can be written as Sj = Sj^ © Wjg © • • • © Wj-i. For any uj G Sj, it follows that 

j-i 

Wj = X] ^J^k4’J,k = ^ cjg^k(l>Jo,k + X] X] (2-21) 


Denote Cj = {cj^k)keAj, dj = {dj^k)kS:Vj and dj = (cj^,dj^,...,dj_i). Then there exists a fast wavelet 
transform (FWT) between the single-scale and the multiscale representations, i.e.. 


cj = Mdj, 


which can be performed with the cost 0(2“^); here 


M = 


Mj_i 0 
0 Ij-i 


Mj-2 0 
0 Ij-2 


Mj, 0 \ 

0 /jo ) 


( 2 . 22 ) 


(2.23) 


and Mj = 


3. Uniform Schemes 


The nonlocal property of fractional operator makes the matrix of its discretizations inevitably dense. We will 
show that the chosen bases being the dilation and translation of one single function render the matrix to have a 
special structure, which greatly reduces the cost of computing and storing the entires. In this sense, these kind 
of bases are superior to the other possible bases, such as the usually used finite element or spectral polynomial 
bases. Meanwhile, based on the benefits of these bases, a simple diagonal preconditioner and the fast transform 
are presented to enhance the effectiveness of the widely used nonlinear or linear iterative schemes. 

We first consider the BVP of dD with (7 = 0. It has the variational formulation: Find u G ^^(D) with 
a=l — /3/2(0 < ,5 < I), such that 

a{u,v) = {f,v) VuGHo(O). (3.1) 

More precisely, using integration by parts and the adjoint property of fractional integral operator leads to 

a{u,v) = (^-Ki 3 D{p oD~^ + {I - p) ^D^^)Du, v'j (3.2) 

= Kp (p qD~^Du + {1 - p)xD^^Du, Dv'^ 

= Kp (p oD-^^^Du, ^D-^''^Dv^ + Kp ((1 - p),D;'^^^Du, . 

The bilinear form a(-, •) : 'Hq (fl) x 'Hq (fl) —>■ K is continuous and coercive [TT], i.e., 

|a(M,'i;)| ~ ||m||c||u||„, a{u,u)^\\u\\l. (3.3) 

For / G ^ 2 ( 0 ), Eq. (13.11) admits a unique solution. Letting Sj be a subspace of Hq(D) with order d, the 
Galerkin approximation uj belonging to Sj satisfies 


a{uj, vj) = (/, vj) V uj G Sj. 


(3.4) 
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If u is sufficiently smooth, following (I2.12|) and the Ced’s lemma, one gets 


inf |1 m-uj||„ - 2 '^(“ ||w||H<i(n)■ 

vj&Sj 


(3.5) 


For space discretization, one can either use the scaling basis or the multiscale basis 


generating the following linear systems, respectively. 


= Fj, 
Ajdj = Fj, 


(3.6) 

(3.7) 


where Aj = a{^j, d)j), Fj = (/, d)j), Aj = a('I'*^, 'I'*^), Fj = (/, dt"^). There are the following Lemmas. 

Lemma 3.1. Let (j){x) G Hq (O), supp^(a;) = [0, d] and 4’J,k{x) := x — fc), 0 < fc < 2*^ — d, /c G N. Then 

a Ajm)=^ {^J,k[ , (t>j,ki, ) if and only if k 2 - h = k '2 - k[. 

Proof. For (p{x) G Hq (11), there holds 

[ [ {x - sf~^ (j)'j,,As) ds dx 

Jo Jo 

-J I J. 

/ {x — s)^~^ (j)' (2'^s — fci) ds (j)' (2’^x — ^ 2 ) dx 

Jo 

(2“*^ (/c 2 + x) — s)^ (j)' (2*^3 — ki) ds (j)'{x) dx 


r(/?) 

23J ^2 {d-\-k2) 

W) 

2'iJ fd |•2~'^{x+k2) 

W) 


'0 ^0 


■)2Ja pd px-\-k2 — ki 


m 


*0 J —k\ 


2‘2.Jcx pd px-\-k2 — ki 


m 


{k 2 + X — s — ki)^ (s) ds (/)'(x) dx 
(x — s + ^2 — ki)^~^(j)' (s) ds (j) {x) dx, 


'0 JO 


which just depends on the value of k 2 — fci. The second part of (13.2|) can be expressed by its first part, i.e.. 




Then the desired result is obtained. 


(3.8) 


□ 


Lemma 3.2. Let(j){x) and (j)j^k{x) be given as above, and 4>{d/2 — x) = 4>{d/2+x). Define 9j^i{x) := 2'^i^9i{2'^x) 
and 9j^i{x) := 2'^/^di(2'^(l — x)) with 9i{x) G 'Hq{TI) and suppdi(x) = [0, d^], where 0 < di < d and i = 1,2. 
Then 

(3.9) 




2J-d-k, . 


Proof. Similar to Lemma l3.11 it follows that 


r(/3) 

22 Jo 

r^) 


/ f {x-sf ^ 9'j As) ds (l)'jj,{x) dx 
Jo Jo 


'0 do 

pd px-\-k 


n X-\-K 

{x k — s)^~^6l (s) ds (j)'{x) dx. 
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By the properties of symmetry and compact support, there exists 




-d-k. xD, 


1 


r(/3) jQ jQ 

22J rO 

Jd. Jo 

Q^Ja pdi nmin{2'^ —k,d} 


m 


{x - sf ^ (j)'j^ 2 J-d-kis) ds 0'J^{x) dx 

(l — s — 2~‘^x)^ ^ 4>' (2*^5 — 2“^ + d + fc) ds di{x) da: 
(s + fc — a:)^ ^ 4>' {s) dsOl (x) dx 


'0 max{0,fc —A;} 


■)2Ja I'd nx-\-k 


pa px-\-K 

/ / (cc + fc — (s) ds (/)'(a;) dx, 

Jo Jo 


r(/?) ao ao 

where the Fubini-Tonelli theorem and min \ — k,d^ = d are used. 

It is also easy to check that for any ii,i 2 G {1,2}, 


□ 

(3.10) 

(3.11) 


Now, from the structure of $ j [55] and the above lemmas, one knows that the matrix Ai := [oDx 
has a quasi-Toeplitz structure, that is, it is a Toeplitz matrix after removing very few rows and columns near 
the boundaries. More precisely, for d = 2, it is a full Toeplitz matrix, but for d = 3 and d = 4, they have the 
following structures, respectively, 


/ 

Oi 


r(a2)^ 

0 

\ 



, 

ai 

H(2-’-2)x(2J-2) 

a2 


5 


V 

02 


r(ai)^ 

Ol 

/ 2-’ 

x2'l 


/ 

Ol 

a2 

r(ai)'^ 


0 

0 

\ 


03 

0-4 

?'(a2)^ 


0 

0 



as 

a4 

H{2-'-0)x{2-’ 

-3) 

a2 

ai 



05 

06 

r{si^Y 


04 

02 


V 

07 

05 

rl&zV 


03 

Ol 

) 


(3.12) 


(3.13) 


where ai are real numbers; ai are vectors, ^(ai) the reverse order of a^; and H^xn is Toeplitz matrix. 

The fact that the bases are obtained by dilating and translating of a single function and the symmetry of the 
bases are essential for obtaining the above results, recalling that they do not hold for the general finite element 
(except linear element) and spectral methods. For the high order hnite difference methods, the similar results 
can be got after modifying the approximation near the boundary for recovering the desired accuracy |5lH5]. 
but it seems that the general theoretical results (stability, convergence and so on) are hard to obtain. Further 
results for the generated matrix are 


=0 Vfc2 - fci < -d. 


(3.14) 

(3.15) 
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Remark 3.1. The structure of <I>j also allows one to compute its Riemann-Liouville fractional derivative 
easily, which can greatly reduce the computational complexity of generating the differential matrix in Galerkin 
and collocation methods (see Section 5). As an example, we present the techniques for d = A in Appendix, 
being similar for other values of d. Then combining with iS.lSl) and i3.15\) . the left differential matrix Ai = 
can be calculated exactly or numerically with the cost 0{N), being superior to the traditional 
finite element and spectral approximation with the cost 0{N‘^) 

3.1. Multilevel Preconditioning 

For an algebraic system with dense matrix, a well convergent iterative method generally has the computational 
cost 0{N \og{N)) or 0{N‘^), which is much less than the cost 0{N^) of the direct method. Moreover, a well 
conditional number and ‘bunching of eigenvalues’ usually bring good numerical stability and fast convergence 
speed mm- In general, for a linear system Ax = 6 , a satisfactory preconditioned system Bx' = b' should have 
the property 


||.B|| < C, ||i? ^11 < C, C is a moderate-sized constant independent of TV; 

and the computational cost for the preconditioning step is cheap. Here, both the matrix Aj and Aj are dense, 
and their condition numbers are of order 0(2^'^“); see Table [7] for Example 5.2. But with the aid of wavelet 
bases, by the norm equivalence, a simple diagonal scaling can lead to a good preconditioned system. In fact, 
define 

V"'= Ajo U Vjo U ■ • ■ U Vj_i, (3.16) 

K = diagf ..., 2"'^““,..., 2-‘^«“,..., 2-('^-i)“,..., 2 -('^-b“ ) . ( 3 . 17 ) 

V'--V-^ '-V-^ ^-V- '' 

#Ajo #Vjo #VJ„1 

Combining the ellipticity (13.31) . norm equivalence (12.131) . and the Riesz representation theorem, one gets that 
for all X £ ^ 2 (V*^), 


\\KAjKx\\i^^^j) = sup 

yehiV-J) 


{kAjKx, y)i^iyj) 


= sup 


\\y\\i2(v-’) 
a{xA'K<^ff y'^K<S>A < \\x'^K^^\\ 


\y^K^‘^ 


|a < 


WVWhiV-’) 


|y|u2(V'') 


X 




\\KAjKx\\i^(^^j-) 


> 


\x^K<i> 




Therefore, there exist Ci,C 2 not depending on J such that 




Ci||a:||/2(v-') < \\KAjKx\\i^^yj) < C'2||a;||i2(v-^)- 


Now, one arrives at 


\\KAjK\\ ^ C 2 , \\{kAjK)-^\\ ^ (1/Ci), 
cond 2 (KAjKuj) = \\kAjK\\\\{kAjK)-^\\ ^ (C 2 /C 1 ). 


(3.18) 


(3.19) 


(3.20) 

(3.21) 


The norm equivalence implies 2 ^A^ 30 one can also define matrix K by the inverse square root 

of the diagonal of Aj, and (13.191) and (13.211) still hold. Usually the current K performs better since it uses the 
information directly from the stiffness matrix, and we will use it in Section 5. Moreover, the cost of generating 
K is only 0{J)-, this is because that by using the translation property of the inner wavelet on the same level. 
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one just needs to calculate the entries a{ipj^k,'4^j,k) near the boundaries and one in the inner part without the 
necessity to assemble Aj. 

Now, one can rewrite (1X71) as the two-sided preconditioned form 


KAjKK-^dj = KFj. 


(3.22) 


Further using (I2.22|l . one gets that 


KM^AjMK R-^M-^cj = KM'^F. (3.23) 

A straightforward product of Aj or Aj to a given vector needs a computational cost 0(2^'^). But if one uses 
the quasi-Toeplitz structure of the matrix, the computational cost can be reduced to 0{J2'^). In fact, one can 
rewrite Aj as 

Aj = diag(A:i)A/ -b diag(A'2)Ar, (3.24) 

where Ki and K 2 denote the coefhcient vectors, formed by the coefficient of space fractional derivative taking 
values at the discretized intervals, and Ai and A^ are quasi-Toeplitz matrices. Using the FFT to the matrix- 
vector product makes the computational cost 0{J2'^) (34]. Finally, because the FWT (having the matrix 
representation M or M^, which denotes the primal reconstruction or the dual decomposition m) can be 
implemented with the cost 0(2'^), if the CG scheme (symmetric) is applied to (13.231) or to the corresponding 
normal equation (asymmetric), the well conditioned number of the matrix implies that the convergence rate is 
independent of the level J; then we can solve it with the total operations 0{J2^). For the general iterative 
schemes, such as GMRES or Bi-CGSTAB, usually one can show that the system with clustered spectrum and 
well conditioned number after preconditioning has an accelerated convergence. What’s more, compared with 
the most existing preconditioners which require the solving of a linear system (see, e.g., the ILU [18] and the 
Strang [15]), the wavelet preconditioning operation reduces to the matrix-vector product, where FWT can be 
used. 

3.2. Multiresolution Multigrid Method 

The multigrid method based on the finite difference discretization for solving the fractional IBVPs have been 
developed in [am], where the transition operators (restriction and prolongation operators) between the grids 
are chosen as the full weight and interpolation operators. In this Subsection, we investigate the MMG method 
for solving fractional IBVPs. We will show that the transition operators in the MRA background can be more 
straightforwardly defined. And using the techniques presented in the above content, the MMG scheme can also 
be fast implemented. 

Denoting Sj as the subspace, Aj : Sj Sj with {AjU!j,Vj) := a(ujj,Vj) \/vj £ Sj and Qj : —?► Sj with 

{QjP^Vj) = {p,Vj) \/vj £ Sj, one arrives at the semidiscrete form; Find uj(t) G Sj,t > 0 such that 


\ Mj(0) = Uj £ Sj. 

Taking the time mesh as 0 = to < < •' • < tv-i < In =T and the stepsizes At„ = —1„, n = 0,..., N —1, 

one gets the backward Euler multiresolution Galerkin method (B-MGM) 

iU]+\vj) + Atr,a{U]+\vj) = iUy + Atr,f{tn+i),vj) WvjGSj; ( 3 . 26 ) 

and the Crank-Nicolson multiresolution Galerkin method (CN-MGM) 

+ uy 



{dU]+\vj) + a 


2 


{f{tn+l/2),Vj) yVjGSj, 


(3.27) 
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where dUj^^ = — C/y)/At„. Introduce the bilinear form := {u,v) + \Atna{u,v), where 

A = 1 for the B-MGM and A = 1/2 for the CN-MGM. Define : Sj —>• Sj with pj,Vj) = Bn+i{pj,Vj) 
Vuj S Sj, and the operator ■ "Hq (D) —>• Sj with Bn+i{P^~^^p,Vj) = Bn+iip,Vj) \/vj £ Sj. Then the MGM 

schemes can be rewritten uniformly as the form 

= g^+\ (3.28) 

where := Uj + AtnQjf{tn+i) for the B-MGM and := -^^AjUj + Uj + AtnQjf{tn+i/ 2 ) for 
the CN-MGM, respectively. Suppose that Uj = ^J,k4’J,k £ Sj, and define cj,gj £ (cj)fc := 

Cj.fe, {gj)k ■■= {gjAj.k), k £ Aj. Denoting = {Bn+i{(j)j,k,<l)j,i))kj^Aj^ algebraic 

representation of (13.281) given by = gj®"®. 

The basic iteration algorithm for the operator equation (13.281) is [55] 

jjn+l.l+l ^ jjn+1,1 ^ J^n+l ^^n+l _ j^n+ljju+l,iy Z = 0 , 1 , 2 ,... (3.29) 

with the error propagation operator /C”®"® ■= I — i?j®’®Sj®'® and the iterator i?”®"® : Sj —> Sj. For the damped 
Richardson and Jacobi methods, the iterators are, respectively, given by 

R^j+^g = ivaiB'J+Y" Y. (9,^J,k)<l^j,k ^g & Sj; (3.30) 

feeAj 

.R"^®5 = ^ ^n+i {(l>j,k,<^j,k) ^ ig,<^J,k) 4>J,k V 5 £ S'j; (3.31) 

fcGAj 

they can also be regarded as the correction with subspaces decomposition Sj = X^fceA^ with = spanjt/j^}; 
and one can also rewrite the Jacobi iterator as 

feeA„- 

where p”+®A- . yk i?„+i(P"®'®’^Aj, (/)j_fc) = Bn+i{vj,4>j,k) ^Vj £ Sj. Now, the MMG D-cycle 

algorithm for (I3.28P reads 

+M”+® (g”+®-S”+®D”+i’') , Z = 0,1,2,... (3.33) 

and the multigrid iterator Alj®"® is defined in Algorithm |T] by induction, where RJ^® and A^"''"® : Sj —>■ Sj are 
defined in the same way as i?"^® and /C"^®. Obviously, the MMG error propagation operator satisfies 

/-M^+i®S"+i® = (/-Af”+®QjS"+/) (3.34) 

"-V-" 

I 

'-V-' 

n 

where the relation QjB'^Xi = S"^®P"®"® has been used; I just denotes the usual two-grid error propagation 
operator; and mi{j + 1) (or m 2 {j + 1)) means that the iterative times mi (or m 2 ) may depend on the level 

j + 1 - 
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Algorithm 1 MMG V-CYCLE ITERATOR 


1 : Fix t = tn+i; for j = Jo, define Assume that 

g G Sj, define the iterator : Sj Sj through the following steps: 

2 : (1) Pre-smoothing: For Xq~^^ = 0 G Sj and 1 = 1,..., mi{j), 




{g - 


3: (2) Coarse grid correction: 


.y.n+1 _ n+1 


+ (5 


-B. 


■n+l^n+1 


”ii 0 ) 


4: (3) Post-smoothing: For I = mi{j) + 2,..., mi{j) + rn 2 {j) -t 1, 

{g - B';+\^+^) 


5: Define 

J a rni{j)+m2{j) + l 


Sj-i —?> Sj-i is defined. For 


Recall that the refinement relation (I2.20p . one can get the prolongation matrix Mj^ straightforward, and it 
holds 

r = M,,oC-+i Vt/;VV = C^r'e5,c5,+i; 

{ - -- _ (3.35) 

where the meaning of is defined after Eq. (13.281) . This means that the transpose of Mj^ is just the 
restriction matrix. Noticing that VC/”^^ G Sj, there holds 




»n+lrrn+l 

j j 




,n+l 


(3.36) 


i.e., 

^n+i ^ mJoB;+>4o, (3.37) 

which actually is the Galerkin identity, facilitating the convergence analysis, but this is not true for the difference 
method. Note that the quasi-Topelitz structure of makes it feasible to be generated directly with the 

cost 0{2^). Using the fast algorithms (FFT and FWT), the matrix-vector product can preformed with the cost 
0{j2^). So the total computational count per MMG step is 0[J2'^) and the storage cost is 0(2*^). 

In the following, we present the convergence analysis of the MMG when A is a Riesz derivative and mi(j) = 
^ 2 {j) = ™o; see Algorithm 1. It is easy to check that B'-^^ is symmetric, and are A-orthogonal 

projectors, and and I — are A-selfadjoint; all of them are considered with respect to P„+i(-, •). 

Lemma 3.3 (see [3]). Assume that : Sj —^ Sj is symmetric with respect to (•,•), positive semi-definite, 

and satisfies 

Bn+i {JC'J~''^Vj,Vj) > 0 Vuj e Sj, 

v„v,) < eP„+i {v„v,) Vu, g{I- p;_Y) S,. 

Then we have 

0 < P„+i ((/ - vj ,v,)< 6Bn+i {vj, V,) Vu, e S,, 


(3.38) 

(3.39) 


where S = e/{e-\- 2mo). 
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Since I — is A-selfadjoint, (I3.39P actually means that its spectral radius 


(/ - = \\I- M^j+^B^+^ I 


= sup 
o^veSj 


Br,+ l{{l-M^j+^B^+^)v,v) 
Bn +1 iv,v) 


< 6 < 1 . 


(3.40) 


When w S [co, 1], 0 < Co < 1, the Richardson method obviously satisfies the requirements of Lemma ESI Since 
the damped Jacobi iteration converges under the condition 0 < w < for any Vj € Sj, there 

exists 


i?„+i {JC]+\j,JC]+^Vj) 

= R„+i {vj,v,) - 2wR„+i {R]+^B:^+\,v,) +a;2R„+i R]+^B^+^Vj) 

= i?„+i {v„v,) - 2u; {R]+^)iB-+\) 

< R„+i {vj,vj) - w (2 - u;a{R;+^B^+^)) [iR]+Y^B^+\, iR]+^)^B^+\) . 

Then it is sufficient to take 0 < w < l/a{R'J'^^B'J~^^) for getting the first condition in (I3.38p . 


feGA,- 


keAi 


— , / 'Y^ Rn+l {£j~^ ^"+1 i^i,k4>j,k,Cj^k4'j,k) 


keAj 


keAj 




feeA,- 


fcGA, 


Y ) ^Uj, Uj) / Bn+1 {Cj,k4'j,k^ Cj,k4'j,k), 

’ y fcGAj 


where v, = E G 5, and £^+^ := P"+'’'=(S"+i)-i 


k(pj,kTCj^k4’3,k) ^ ^ 'Y^ ll'^t,fc'/'i,fc|lL2(n) ^ lkillL2(n)’ 


(3.41) 


fceAj 

By Bernstein estimate and uniform stability given in (12.101) and (12.81) . one gets 

feGAj 

By the Aubin-Nitscale trick [Sl fTUlfT^ . there exists 

|(/ - PYiH^L^in) - 2-^^“ ||(/ - ((/ - PYi>3^ (I - PY>^) ■ (3-42) 


Y^ Pn+l {Cj, 

fcGA,- 


IW“(a) 


Noting that (/ — P^^i)vj G Sj, then the second requirement in p.38l) holds. The proof of the convergence for 
MMG with Jacobi iterator is completed. 

4. Multiscale Adaptive Schemes 


Although there are works to discuss the low regularity (especially the weaker regularity at the area close to 
boundary) of the solutions for fractional PDEs, it seems few of them are for designing the adaptive algorithms, 
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which have been well developed for classical PDEs. In the classical finite elements approximation, adaptivity is 
usually driven by so-called local a posteriori error estimate (an efficient and reliable error indicator consisting of 
local terms and being easy to compute). In the following we first show the challenge when using the traditional 
finite element method to adaptively solve the fractional BVPs, then numerically demonstrate the performance 
of the wavelet adaptation algorithm. 

For the linear element and /3 £ (1/2,1), with a given positive integer K, we define the mesh 

0 = xo < a;i < • --xk-i < xk = h = (xi_i,Xi). 


Then the posteriori error for the BVP is bounded by 

K 

a{u -Uh,u- Uh) ~ y] hj°‘ \\f + Kji {pQDl°‘uh + (1 - p)xDl“uh) , (4-1) 

where hi = xt — Xi_i, and /3 £ (1/2,1) ensures that oDl°‘uh and xDl°‘uh belong to L 2 . 

Now we prove mi- Denote by the operator for the piecewise linear interpolation associated with {xi} 
and eh = u — Uh- For any v £ there exists 

a{eh,v) = a{eh,v - 11,*^) = (/,u - n,ju) - a {uh,v - n,ju). 

Combining (13.2L {v — Whv){xi-i) = {v — n/iu)(xi) = 0, and the regularity of Uh leads to 


',{eh,v) = T [ {f + k/ 3 {poDl°‘uh +jl - p)xDl°‘uh)) {v-Hhv) dx 

K 

^ ^ 11-^ + + (1 - p)xDl°‘uh) 11^^ - nhu|| 

K 

^ 11/ + '^^ {PoDl°‘uh + {1 - lkllw-(/,) 


K 


\ 11-^ + {poDl°‘Uh + (1 - p)xDl°‘uh)\\l^^j.) |k||«c(n). 

\ i=i 


Then (14.111 follows by taking v = Ch in the above inequality and using a(?;, u)'^||u||^c«(q). The term || ■ \\h 2 {ii) 
involves nonlocal calculations; so it can not be used directly as a local error indicator. The following in Algorithm 
[21 we will show that for the wavelet methods of fractional BVPs, the local regularity indicator can be the wavelet 
coefficient when the (to be determined) solution is represented by the multiscale bases: small coefficient implies 
good local regularity while big one indicates the opposite. One of the main features of Algorithm [2] is that the 
finest grid resolution can be automatically determined by the given tolerance e(j). In order to gain a more robust 
and faster algebra solver, has been scaled by the inverse square root of a {4’j,k,'4’j,k), and the multiscale 
approximation of the solution at the current scale has been as an initial guess for the iteration in the finer 
scale obtained after adding the wavelets. For constructing the refined index set, thanks to the tree structure 
of wavelet singularity detection, we first include a coarsening step by thresholding the latest available wavelet 
coefficients to get a significant index set, then add all their children. If j = Z -|- 1 and k £ {2A, 2A -I- I}, then 
the wavelet indexed by (j, k) is called a child of the wavelet indexed by (Z, A). One can further extend the index 
set by including the horizontal neighbors of the wavelet indices already included. Such an extended index set 
associated with the index (Z, A) is called an adjacent zone, which is denoted by A/,a- In Algorithm |21 the index 
set is continuously updated to resolve the local structures that appear in the solution. One can dynamically 
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Algorithm 2 ADAPTIVE WAVELET SOLVER FOR THE BVP 

1 : Given e{j),Itmax, Jo 
2: m = 0 

3: Solve the equation in space Vj^+i to get the initial approximation coefficients (cj , dj ) and the index 
A- = (Jo,A),Ae 


4: 

5: 

6 : 

7: 


9 

10 

11 


12 

13 

14 


repeat 

Determine the significant index set A by e(j) 

Check the adjacent zone index set A/i,A of each (I, A) € A; denote Ajy = U(i^A)GAA/i,A and establish 
Am+i = A U Ajy 
for (I, A) G do 

= / a,A)€A- 

1 otherwise 

end for 

Jo Jo 

Solve the algebraic matrix equation, resulted from the discretization in the nonlinear approxima¬ 
tion space Vjg+m+i(^) C Vjp+m+i(0), by appropriate iterative scheme with the initial guess 

/* m+l * jm+1 * jm+1 '> 

V ^Jo ’ ’ ■ • ■ ’ ‘^Jo+m+ll 

Determine A = {{I, A) : (l, A) G A™+^, > e(/)} 

m = m -b 1 

until m > Umax or A = $ (empty set) 


adjust the number and locations of the wavelets used in the wavelet expansion, reducing significantly the cost 
of the scheme while providing enough resolution in the regions where the solution varies significantly. The m-th 
approximation of the solution is given by 


Ujg+m = ^ ( 4 . 2 ) 

(J,A)GAjgUA™ 

where Ajg U A™ is the irregular index set, and represents the normalization of t/^i^x- Finally, after getting 
the sufficiently accurate approximation, the corresponding single scaling representation can be got by the FWT. 

The following we develop adaptive algorithm for time-dependent problem. With the time partition 0 = to < 
ti < ■ ■ ■ < tN = T and stepsizes At„ = tn+i — tn, n = 0,..., A — 1, for some time r G tn], one arrives at 
a system of the form 

dtu\r + Au{-,t) = f{-,T). ( 4 . 3 ) 

The numerical solution [/" can be uniquely represented (a unique decomposition) in one of the subspaces of 
Sj :=SjgUWjgU---Wj_,: 


j-i 




keAjg 


which is equivalent to the unique coefficient vector 


j=jQ AeVjn^" 


P” := (E}o>dX,...,d:_i). 


( 4 . 4 ) 


(4.5) 


For establishing the algebraic system of p", one can still use the Galerkin scheme; as an extension, the collocation 
method based on the semi-interpolation wavelet (see eq. (12.16I) - (I2.191) 1 can also be considered, which replaces the 
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test function a} U {V'j.a} , (j, A) G of the Galerkin scheme by the Dirac distribution S centered at Xi, be- 
ing the collocation point corresponding to the index set f?", a subset of U {(2^ + k=o j=Jo' 

Generally, the collocation method is more convenient and efficient for problems with variable coefficients and/or 
nonlinear terms. Details of such a scheme are provided in Algorithm |31 the steps are very similar to Algorithm 


Algorithm 3 ADAPTIVE WAVELET SOLVER FOR THE IBVP 
1 : Given time partition {t„}, Jmax and threshold e(j) 

2 : Gonstruct the initial irregular index set and the multiscale coefficients by uq 
3: for n = 1,..., A do 

4: Based on to solve p” on the index set 

5: Threshold p" to obtain the significant index set 

a” := {(J/A) e : |4 a| > e{j)} 

6 : Add the adjacent indices to it, and denote the result by 

7: if and are different then 

8 : For every index (j. A) S C/” not in the corresponding wavelet coefficients are initialized with 0, 

and denote the result by p” 

9: end if 

10: end for 


[21 except that for treating the structures appearing in the solutions as they evolve, the computational index 
needs to dynamically adapt to the local change of the regularity of the solution. For an implicit or explicit 
time integration, we use wavelet amplitudes of the approximate solution at the current time level to construct 
the irregular index for the approximate solution of the next time level. The initial irregular index can be 
constructed by adding the adjacent zone to the significant index set of the initial solution u{x,0) = g{x). 


5. Numerical results 


In order to illustrate the accuracy and efficiency of the proposed numerical schemes, we apply them to solve 
the BVP and/or IBVP (II.II) . Example 15.11 is used to discuss the implementations of the MGM for the BVP 
and the collocation method for the IBVP, and in particular the convergence orders are carefully verified. We 
use Example 15.21 to show the powerfulness of the provided multilevel preconditioner and MMG. And Example 
15.31 is used to illustrate the effectiveness of the presented wavelet adaptive schemes. 

Example 5.1. 

Gonsider the MGM for the BVP (II.II) with <7 = 0 and p = K/j = I, and the source term 


f{x) 


2a^^ _ r(u+I) 

r(/3 + I) + 


The exact solution of the problem is u{x) = x'^ — x^ . It is well known that if u > 0 and u ^ N, then 
u G For /? = 4/5, the numerical results are listed in Tables |T] and |21 which confirm that if the 

analytical solution is smooth enough, the convergence order is d and d — a in the L 2 and ?^“-norm, respectively. 
Otherwise the convergence order is limited by the regularity of the solution, but the approximation accuracy 
is improved when the high order bases are used. Moreover, If the modified Galerkin method, e.g., the one 
proposed in [14], is used, for this type problem one can have a convergence rate d — /3 for the sufficient smooth 
source term /; when / = 1, the numerical results are listed in Table O We want to emphasize that including 
the boundary bases is very important to ensure the polynomial exactness (known as the Strang-fix condition), 
which is the foundation to have the desired convergence results. The numerical results in Table [4] are for the 
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Table 1. Numerical results of the BVP dm, solved by MGM, with q = 0^ p = np = 1, and /3 = 4/5. 


d J 

1/^4 

V = 17/10 


L2-Err 

L2-Rate 

a{u — uj, u — 

7/“-Rate 

L2-Err 

L 2 -Rate 

6 

17589e-04 

— 

8.5677e-04 

— 

1.4535e-05 

— 

d^2 7 

4.3968e-05 

2.0001 

3.1635e-04 

1.4374 

3.6314e-06 

2.0009 

8 

1.0993e-05 

1.9999 

1.1829e-04 

1.4192 

9.0287e-07 

2.0079 

6 

6.2317e-07 

— 

6.2807e-06 

— 

1.0342e-06 

— 

CO 

II 

7.7779e-08 

3.0021 

1.1896e-06 

2.4004 

2.2509e-07 

2.2000 

8 

9.7152e-09 

3.0011 

2.2531e-07 

2.4005 

4.8988e-08 

2.2000 


Table 2. Numerical results of the BVP dm, solved by MGM, with q = 0^ p = Kp = I, and [3 = 4/5. 


J 

V = 11/10 

u = 21/10 


d = 3 

d = 4 

d = 3 

d = 4 


Z/2-Err L 2 -Rate 

Z/2-Err L2-Rate 

L2-Err L 2 -Rate 

Z/2-Err Z/2-Rate 

6 

1.4385e-05 — 

8.0390e-06 — 

1.2656e-07 — 

3.2703e-08 — 

7 

4.7453e-06 1.6000 

2.6516e-06 1.6002 

2.0865e-08 2.6007 

5.3930e-09 2.6002 

8 

1.5654e-06 1.6000 

8.7469e-07 1.6002 

3.4407e-09 2.6003 

8.8950e-10 2.6000 


Table 3. Numerical results of the BVP (II.IL solved by MGM, with g = 0, p = k /3 = 1, / = 1, 
d = 4, and p = 4. 


J 

p = 

4/5 

p = 

1/2 

/3 = 

1/5 


L2-Err 

L2-Rate 

L2-Err 

I/ 2 -Rate 

L2-Err 

L2-Rate 

6 

2.4209e-07 

— 

9.0406e-09 

— 

3.0206e-10 

— 

7 

2.6360e-08 

3.1991 

8.0061e-10 

3.4973 

2.1733e-ll 

3.7969 

8 

2.8694e-09 

3.1996 

7.0772e-ll 

3.4999 

1.5568e-12 

3.8032 


cases that the boundary base functions are absent (one base is removed for d = 3 and two for d = 4). At this 
moment the exact solution after zero extension is required to have sufficient regularity to recover the desired 
convergence order. The similar observations are also detected for the finite difference methods, and the ways of 
recovering the optimal convergence orders are presented in [51l42j. 

Table 4. Numerical results of the BVP dni), solved by inner MGM, with q = 0, p = np = 1, 
and /3 = 4/5. 


d 

J 

u{x) — 

1 

to 

u(x) — x^ 

(x-lf 

u(x) — x^ 

G-i)" 



L2-Err 

L2-Rate 

L2-Err 

L2-Rate 

L2-Err 

L2-Rate 


6 

9.8488e-3 

— 

1.2187e-06 

— 

5.9382e-07 

— 

3 

7 

4.9332e-3 

0.9974 

1.5229e-07 

3.0004 

7.5027e-08 

2.9845 


8 

2.4688e-3 

0.9987 

1.9027e-08 

3.0007 

9.4230e-09 

2.9931 


6 

1.9092e-2 

— 

8.9636e-05 

— 

5.9852e-08 

— 

4 

7 

9.6958e-3 

0.9775 

2.2395e-05 

2.0009 

3.7671e-09 

3.9898 


8 

4.8857e-3 

0.9888 

5.5941e-06 

2.0012 

2.3616e-10 

3.9956 


We further consider the collocation method for the variable-coefficient version of the IB VP (11.11) . The 
collocation points are chosen as {l/2'^+^, fc/2'^|^_j^^, 1 — 1/2'^+^}, and the approximation properties of the cubic 
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spline collocation method are discussed in the space Sj := span{<i>j} with <i>j = {(j)j^k,k G Aj,ci = 4}. The 
considered equation is 


Ut-{hx^-%Dl-Pu + k2{l-x)^-^ ,dI-^u) = f te{0,T] (5.1) 

with the right-hand term 

f{x,t) = -12exp(-t)|x^(l - -I-i -I-^2(1 - 

+ (J3 + 1K„ + 2) + '•'^<1 - “^>■‘1 } 

and the initial condition m(x, 0) = x^{l — x)^. Then it can be checked that the analytical solution is u{x,t) = 
exp(—t)a;^(l — xY. 

The Crank-Nicolson scheme is used to get the full discretization approximation of with the time stepsize 
1/2^'^ and T = 1/2. Table O shows the expected convergence order 2 + (3 oi collocation method for fractional 

Table 5. Convergence performance of the cubic spline collocation method with fci = ^2 = 1- 


J 

/3 = 

2/10 

0 = 

8/10 

0 = 

0 


Loo-Err 


Lcx3-Err 

Loo-Rate 

Loo-Err 

Loo-Rate 

5 

5.8773e-05 

— 

1.8431e-06 

— 

1.6167e-04 

— 

6 

1.2S62e-05 

2.1920 

2.6580e-07 

2.7937 

4.0533e-05 

1.9958 

7 

2.8036e-06 

2.1977 

3.8223e-08 

2.7978 

1.0441e-05 

1.9990 


PDE, agreeing with the classical conclusion when ,5 = 0. Though the superconvergence can be obtained for 
classical PDE by carefully averaging the derivative values gotten in collocation points, it seems that this result 
maynot be directly extended to fractional PDE. Eor convenience, let’s consider the frequently used Hermite 
spline collocation method. Take the collocation space 


Vj := span {712(2a;-b l)|Q,7ri(2'^a; - fc), 712(2'^a: - fc), 712(2*^3; - 2^^ -b 1)) , 

where denotes the restriction inD, k = 0,1,-- - ,2*^ — 2, and tti, 712 are the cubic Hermite compactly supported 
functions given as 


7ri(a;) = 
712(3;) = 


—a;^(2a; — 3) 

(a; — 2)^(2a; — 1) 

a;^(a; — 1) 

{x — 2Yix — 1) 


0 < a; < 1, 

1 < a; < 2, 

0 < a; < 1, 

1 < a; < 2. 


To determine the unknown coefficients, one need total 2*^+^ collocation points. As well known for classic 
PDE (/3 = 0), when the general collocation points such as the third-quarter points of every interval (* + 

l)/2‘^], 7 = 0,1, • • • , 2"' — 1 are used, the convergence order is 2. But if the Gauss nodes are used, one arrives 
at the superconvergence result of order 4. Unfortunately, in both case the convergence order are 2 -b /3 for the 
fractional PDE, except the approximation accuracy maybe improved. The numerical results are presented in 
Table m where the abbr ‘Equi’ and ‘Gauss’ denote the two types of collocation points mentioned above. 

Remark 5.1. Unlike the Galerkin method, the differential matrix of right derivative in eolloeation method is 
not the transpose of its left twin. Instead, eomhining the symmetry of the sealing spline bases and collocation 
points, it is easy to prove that 


Ar = A;(end : —1 : l,end : —1 : 1). 
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Table 6 . Convergence performance of the cubic Hermite collocation method with fci = 1 and 
k2 = 0 . 


J 

/3 = 5/10, Equi 

/3 = 5/10, Gauss 

0 = 0, 

Equi 

0 = 0, Gauss 


^oo-Err Loo-Eate 

Loo-Err Loo-Ea-te 

Loo-Err 

Loo-Eate 

Loo-Err Loo-Eate 

5 

8.2952e-06 — 

1.9362e-07 — 

3.5693e-05 

— 

3.5875e-08 — 

6 

1.4663e-06 2.5001 

3.2737e-08 2.5642 

8.9270e-06 

1.9994 

2.2506e-09 3.9946 

7 

2.5912e-07 2.5005 

5.6891e-09 2.5247 

2.2316e-06 

2.0001 

1.4098e-10 3.9968 


Example 5.2. 

Now, we focus on the wavelet multilevel schemes for solving the fractional PDEs. The presented numerical 
results with d = 2, and in this case the coefficient matrix has a full Toeplitz structure. The matrix-vector 
product is performed by FFT. For the other bases, the computational procedure is almost the same after a 
slight modification, e.g., when d = 3, AiCj := Cj can be decomposed into several 

blocks with H being the Toeplitz matrix: 

{AiCj){l)= [ai,r(a2)^,0]cj, 

{AiCj) (2 : end — 1) = Cj(l) ai + Hcj{2 : end — 1) -b Cj(end) a2, 

(^iCj)(end) = [a2,r(ai)^,ai]cj. 

For the BVP, we first reveal that the multilevel preconditioning brings a uniform matrix condition number 
and an improved spectral distribution. Considering the BVP (11.11) with Kp = l,p = 1 and Kp = l,p = 1/2, 
the condition numbers for different /3 are presented in Table [TJ one can see that without preconditioning, 
the condition number of the stiffness matrix behaves like which means the conditional number 

increases fast with the refinement especially when j3 is small. After preconditioning, the uniformly bounded 

Table 7. Primal condition numbers of the BVP CH) with 9 = 0, = 1, and d = 2. 


J 

p = 1, /3 = 

1/2 

p= 1/2,/3 

= 1/2 

p = 1 . /3 = 

1/5 

p= l/2,/3 

= 1/5 


Con-Num 

Rate 

Con-Num 

Rate 

Con-Num 

Rate 

Con-Num 

Rate 

8 

1.4763e+03 

— 

1.8304e+03 

— 

8.7494e+03 

— 

9.1119e+03 

— 

9 

4.1754e+03 

1.5000 

5.1784e+03 

1.5003 

3.0467e+04 

1.8000 

3.1732e+04 

1.8001 

10 

1.1810e+04 

1.5000 

1.4648e+04 

1.5002 

1.0609e+05 

1.8000 

1.1050e+05 

1.8000 


condition numbers with different wavelet preconditioners are obtained; see Table [3 where ‘inte-’, ‘Semi-’, and 
‘Bior- (d)’ denote the interpolation wavelet, semiorthogonal wavelet, and biorthogonal wavelet respectively, 
having been introduced in Section 2. Note that when performing the decomposition by semiorthogonal and 
biorthogonal wavelets, the interpolation wavelet has been used for and 82 - We also display the matrix 
eigenvalue distribution for /3 = 1/5 in Figures [T] and [2] they show the preconditioning benefits of a more 
concentrated eigenvalue distribution. 

To explore the effectiveness of this preconditioned system, we numerically solve BVP (HID with 


/ 


1 

r(i + ^) 


{-2x^ + I3x^-^) , 


/ 


1 

2r(i-b/3) 


and 


{-2x^ + - 2(1 - xf + P{1 - xf-^) 
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Table 8. Preconditioned condition numbers of the BVP (II.ID with g = 0, K/j = 1, d = 2, and 

Jo = 0. 


p 

J 

II 

/3 = 1/5 

Inte- 

Semi- 

Bior- (d = 2) 

Inte- 

Semi- 

Bior- (d = 2) 


8 

3.0970 

5.8363 

13.3957 

1.5953 

10.2897 

12.2315 

p ^ 1 

9 

3.2286 

6.1158 

14.4784 

1.6269 

10.5561 

12.9767 


10 

3.3457 

6.3622 

15.4103 

1.6540 

10.7779 

13.5788 


8 

3.2614 

8.0344 

12.7702 

1.5935 

11.1094 

12.2830 

II 

9 

3.4745 

8.1648 

13.6624 

1.6296 

11.4094 

13.0063 


10 

3.6686 

8.2634 

14.4026 

1.6608 

11.6511 

13.5854 



X 10 


X 10 


X 10 


xIO 


Figure 1. Eigenvalue distribution of the BVP matrix with p = 1 (first two) and p = 1/2 (last two). 



Figure 2. Eigenvalue distribution of the preconditioned BVP systems with the interpolation 
wavelet (hrst line) and the semiorthogonal wavelet (second line), respectively (the hrst two 
columns are for p = 1 and the last two columns for p = 1/2). 


for p = 1 and p = 1/2, respectively. We use GMRES and Bi-CGSTAB to solve the algebraic system before and 
after preconditioning, and the numerical results are given in Tables [9l and 1 1 0[ respectively. The comparisons for 
the two methods are made almost with the same L 2 approximation error, not listed in the tables. The stopping 
criterion for solving the linear systems is 


<k)\\i 


< le - 8, 


with r{k) being the residual vector of linear systems after k iterations. It should be noted that the GMRES 
method for p = 1 without preconditioning stops before reaching this criterion. In fact, by the two-dimension 
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FWT and the properties of tensor product, these proposed preconditioner can also easily apply to two dimension 
and it also works well for the algebraic systems generated by the finite difference methods, e.g., [201130] . 

Table 9. Numerical results of the BVP (II.IL solved by GMRES and Bi-CGSTAB, with q = 

0, K/3 = 1, /3 = 1/5, and d = 2. 


J 

p = 1, GMRES 

P = 1/2. 

GMRES 

P ^ 1, Bi-CGSTAB 

p = 1/2, Bi 

-CGSTAB 


Iter 

CPU(s) 

Iter 

CPU(s) 

Iter 

CPU(s) 

Iter 

CPU(s) 

8 

2.5500e+02 

0.3443 

1.1800e+02 

0.0915 

2.6350e+02 

0.0672 

1.1700e-|-02 

0.0303 

9 

5.1100e+02 

1.5217 

2.2000e+02 

0.3230 

5.3550e+02 

0.2408 

2.0950e+02 

0.1012 

10 

1.0230e+03 

7.4783 

4.1200e+02 

1.3219 

1.1665e+03 

0.6731 

3.9150e+02 

0.2280 


Table 10. Numerical results of the BVP solved by the preconditioned GMRES and 

Bi-CGSTAB, with q = 0, Kp = 1, fd = 1/5, d = 2, and Jq = 0. 


p J 

GMRES, Inte- 

GMRES, Semi- 

Bi-GGSTAB, Inte- 

Bi-CGSTAB, Semi- 


Iter 

GPU(s) 

Iter 

CPU(s) 

Inter 

CPU(s) 

Iter 

CPU(s) 

8 

13.0 

0.0094 

27.0 

0.0203 

8.0 

0.0077 

19.0 

0.0198 

p ^ 1 9 

13.0 

0.0115 

28.0 

0.0258 

9.5 

0.0133 

20.0 

0.0272 

10 

13.0 

0.0209 

28.0 

0.0452 

9.5 

0.0149 

22.0 

0.0376 

8 

9.0 

0.0075 

25.0 

0.0227 

6.5 

0.0064 

18.0 

0.0201 

II 

9.0 

0.0084 

26.0 

0.0248 

7.5 

0.0087 

20.0 

0.0267 

10 

9.0 

0.0163 

26.0 

0.0370 

8.0 

0.0126 

21.0 

0.0363 


Secondly, we use the MMG to solve the fractional IBVPs (El) with the exact solution u{x,t) = exp(—— 
x"^), q = Kp = 1, and the suitable source term and initial condition. It can be noted that because of the constant 
diagonal elements of the stiffness matrix, the Richardson and the Jacobi iterations used in MMG are actually 
equivalent. Eor mi(j) = m2{j) = 1, Jq = 3, the numerical results of CN-MMG are given in Tables [TT] and [121 
where ‘Iter’ denotes the average iteration times and ‘CPU(s)’ the computation time also including the time of 
the calculation of coefficient matrix ,j = Jq, - ■ ■ J and the right term. The initial iteration vector at 
is chosen as the approximation at and the stopping criterion is 

< 2 -^/^ X (le-9), 

OO 

where is the approximation vector after the / iteration. Of course, the EET and the EWT are used to 

accelerate the process. ‘Gauss(s)’ denotes the computation time of the Gaussian elimination method; for fair 
comparison, the EET is also used to get the matrix-vector product appeared in the right-hand term at the time 

tn-t-l. 

Eorp = 1/2, the coefficient matrix is symmetric. And if we choose a; < l/Amax,Amax = (o' , 

then the MMG is convergent and the average iteration number is slightly affected by the choice of w. It also 
seems that the restriction to lo can be relaxed to some extent in real computation. When p ^ 1/2, even though 
there are no strict theoretical prediction, the numerical results show when w > Amax, the iteration may be 
divergent; see Table [T^ Here, we get the value of Amax by the Matlab function eigs; it can also be estimated 
by the Gerschgorin Theorem or the Power method. 

Example 5.3. 

In this example, we focus on the previously proposed ad-hoc wavelet adaptive algorithms for the fractional 
PDEs. The BVP is solved by the biorthogonal wavelet bases produced by {d = 3,d = 3), and the IBVP 


n+l,Z n+l,Z —1 

Cj -Cj 
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Table 11. Numerical results of the IB VP (11.11) . solved by CN-MMG, with q = np = 1, p = 1/2^ 
u = 1 , T = 1, and At = 1/2"^. 


UJ J 


13 = 

7/10 



13 = 

2/10 


L2-Err 

Iter 

CPU(s) 

Gauss(s) 

L2-Err 

iter 

CPU(s) 

Gauss(s) 


8 

5.6512e-07 

5.03 

1.0782 

0.2541 

7.7427e-07 

7.97 

1.6402 

0.2692 

4/(SAmax) 

9 

1.3673e-07 

5.00 

3.2092 

1.6474 

1.8554e-07 

7.00 

4.3035 

1.6994 


10 

3.3486e-08 

4.32 

9.2786 

18.5861 

4.1703e-08 

6.03 

12.2416 

18.9016 


8 

5.6512e-07 

4.00 

0.8895 

0.2528 

7.7431e-07 

6.00 

1.2449 

0.2493 

C^Ajnax) 

9 

1.3673e-07 

4.00 

2.6654 

1.6631 

1.8550e-07 

5.01 

3.2196 

1.7100 


10 

3.3488e-08 

3.59 

8.0329 

18.7333 

4.1709e-08 

4.90 

10.3390 

18.4800 


Table 12. Numerical results of the IB VP (ED, solved by CN-MMG, with q = p = up = I, 
P = 7/10, T = 1, and At = 1/2-^. 



J 

uj — 2/(5An;iax) 


^ = 4/(5A 

max ) 

oj — 6/(5Ainax) 



L2-Err 

Iter 

CPU(s) 

Gauss(s) 

Iter 

CPU(s) 

Gauss(s) 



8 

1.2500e-06 

14.79 

2.9692 

0.3550 

10.95 

2.2017 

0.3667 

nc 

vge. 


9 

3.1242e-07 

12.95 

7.6719 

3.0322 

9.80 

5.8060 

2.9866 

nc 

vge. 


10 

7.9268e-08 

10.80 

20.3230 

31.3422 

8.13 

15.4748 

31.0381 

no c 

vge. 


8 

1.7059e-06 

13.47 

2.7015 

0.4289 

10.44 

2.1180 

0.4300 

no c 

vge. 

11/10 

9 

5.0960e-07 

11.21 

6.8185 

3.1898 

9.08 

5.5673 

3.1589 

„„ 

vge. 


10 

1.5759e-07 

9.01 

17.7328 

31.7698 

7.36 

14.8017 

31.9835 

no 

vge. 


by the semi-interpolation wavelet bases. We first consider the BVP (11.11) . the regularity of its exact solution, 
u{x) = (1 — “ (1 ~ 2 :), is weak at the area close to the right boundary; and the parameters = l,p = 0 , 

and the source term 

r( 21 / 10 )(l-x)^-®/i° 

“ r (/3 + i/io) ^ f(^) ■ 

In the algorithm, we take Jo = 3,e(j) = le — 5. For every iteration step, the finally extended irregular indexes 
are obtained by firstly adding the children of all the significant indexes and then including two neighbors, 
i.e., the right and left neighbors, of each index of the extended irregular indexes. When P = 1/2, the sets of 
wavelet indices that corresponding to the adaptively chosen wavelets and the corresponding error u — uj^+rn are 
presented in Figures |3l where the blue bar denotes that we have used all the scaling bases in the coarest level 
Jq. One can see that the algorithm in fact automatically recognizes the whereabouts of the boundary layer of 
the solution u, and adds wavelets locally to there. It also reveals that the newly added computational costs are 
spended in the most needed place, and the large peaks of the errors are successively reduced. Moreover, for 
different P, from the decreasing of the L 2 approximation error of the adaptive and uniform Galerkin schemes 
with the increasing of the freedom N (the loglog coordinate) in Figure IH one can see that the adaptive MGM 
is remarkably superior to the uniform MGM. 

Secondly, we consider the IBVP (11.11) with np = p = 1, P = 5/10, the initial condition u{x, 0) = x'^{l — x), 
and the source term 

f{x, t) = exp(3t)a;2°‘+^(l - a:)(3 4- 20 Inx) -f exp(3t) p( 20 ^+ 3 +^^) 

Its exact solution is it = exp(3t)x^°‘+'^(l — x), which has a strong gradient at somewhere as shown in Figured 
(left). In the computation, based on Algorithm [3] the semi-interpolation adaptive wavelet collocation method 
is used; and the time step At = 1 / 2 ^'^”*“®, Jq = 3,Jmax = 10,e(j) = le — 5, and T = 1. In this adaptive 
algorithm, for every time step, the index extension techniques being used are the same as the ones for the BVP, 


X 


.20t+2+P 


20 t + 5 
20t + 3 + P 


X — 1 
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1st iteration, N=32 






.-eDifference to exact solution 




0.2 0.4 0.6 0.8 


Figure 3. Distribution of adaptive wavelet bases and curve of the approximation error gotten 
by Algorithm [2j 


and the chosen collocation points are just the ones corresponding to the reserved wavelet bases. For t = 1, the 
adaptive solution and the distribution of semi-interpolation wavelets are displayed in Figure [5] Further seeing 
the global picture, Figure [5] (right), one can easily notice that the high level wavelets and collocation points 
mainly concentrate on the area with steep gradient, being exactly as what we have desired. 

6 . Conclusion and Discussion 

This paper focuses on digging out the potential benefits, providing the techniques, and performing the 
theoretical analysis and extensive numerical experiments in solving the fractional PDFs by wavelet numerical 













































































MULTILEVEL AND MULTISCALE SCHEMES FOR FRACTIONAL PDE 


23 




Figure 4. L 2 errors versus freedom N for the adaptive and the uniform Galerkin approxima¬ 
tions with /3 = 1/2 (left) and /? = 4/5 (right), respectively. 


t= Is, Error= 6.6469e-05 t=ls, N=89 




Figure 5. The adaptive solution (t = 1) and the corresponding distribution of the semi¬ 
interpolation wavelets gotten by Algorithm [31 




Figure 6 . Evolution of the solutions and corresponding distribution of the collocation points 
gotten by Algorithm |3l 


methods. The multiscale (wavelet) bases show their strong advantages in treating the fractional operators 
which essentially arise from the multiscale problem. Even the scaling bases also display their powerfulness in 
saving computational cost when generating stiffness matrix, i.e., by using the scaling bases, the stiffness matrix 
has the Toeplitz structure. The way of generating effective preconditioner is presented for time-independent 
problem and multigrid scheme for time dependent problem is detailedly discussed. We numerically show that 
the heuristic wavelet adaptive scheme works very well for fractional PDEs; in particular, it is still easy to get 
the local regularity indicator even for the fractional (nonlocal) problem; and the algorithm descriptions are 
provided. 
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After finishing this work, one of the directions of our further research appears, i.e., applying the wavelet 
compression property to fractional operator. A key difference between the fractional and classical operators is 
that the former is non-local, and then both the matrixes generated by the scaling and the multiscale bases are 
no longer sparse. Fortunately, the wavelet compression not only allows one to obtain a sparse representation of 
functions, but it seems also effective for the fractional operators. Considering the discretization of the operator: 

Au = -D Du 

in the approximation space Sj with J = 10, we first compute the matrix Aj or Aj (here the multiscale wavelet 
bases also have been normalized with D, proposed in Section 3 ). Then we get the compressed matrix by setting 
all entries of Aj or Aj with modulus less than e = 10“^ x 2“'^ to zeros. The comparison results are displayed 
in Table IT^ and Figure[7], where (• %) denotes the percentage of the non-zero entries of the compressed matrix. 
It can be seen that many entries in Aj are so small that they can be omitted to retrieve the famous finger 
structure, whereas essentially all entries in Aj are significant. In the future, we will investigate the effective 
ways of using wavelet compression to get the paralleled sparse approximate inverse (SPAI) preconditioner and 
to perform the low-cost multiscale matrix-vector product. 




Figure 7. Sparsity patterns of the compression matrix Aj by using the semiorthogonal with 
(3 = 5/10 and d = 2 (left) and the biorthogonal wavelet with d = 2 and d = A bases (right). 


Table 13. Compression capacity of the different bases for operator A. 


■0 

0 = 

8/10 

^ ^ 5/10 

p = 

2/10 

Note 

Aj 

Aj 

Aj 

Aj 

Aj 

Aj 

wavelet 

compression 

d — 

2, Inte — 

99.80% 

99.07% 

99.80% 

80.38 

99.80% 

48.04% 

interpolation 

No/Yes 

d - 

2, Semi — 

99.80% 

6.66% 

99.80% 

6.14% 

99.80% 

5.31% 

semiorthogonal 

Yes 

d - 

to 

II 

99.80% 

8.15% 

99.80% 

7.89% 

99.80% 

7.15% 

biorthogonal 

Yes 

d = 

3, d - 3 

100% 

10.18% 

100% 

9.52% 

100% 

8.30% 

biorthogonal 

Yes 


Appendix 

Here we present the techniques for computing the left fractional derivative of the base function = {(/yfc, k £ 
Aj,d = 4}. Besides (/(x) and given by (12.161) and (12.171) . define 

(j)a{x) = 3x+ - ^xl + ^xl - 2{x - l)i^ -f i(x - 2);^. 
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Then x — k)\^_^ , (t)b (2'^(1 — a;)) ,4>a (2'^(1 — x)) } is a Riesz bases of Sj. For 

xo > 0 , it is easy to check that 

oD^r^ (iJ(x - xo)w(x)) = H{x- xo)xoDI~^v{x), 
where H (x) denotes the Heaviside function. By the well-known formulae 

T{u + l){x-aY+P-^ 


aDl-P{x-ar = 


r(u + ^) 


u e N, 


(6 - ax)% = {h- axf -f (-l)''-^(ax - b)%, k S N+, 


for b/a > 0, k G N+, there exist 

oDl-^{ax - b)l = 

k 

oDl~’^{b - axY+ = {-lY-^oDl~^{ax - 6 )+ + ^ 

771—0 

Define 

M,(i) := (o + 1)4 - 34+' + 

+sw3)H<*-‘)+"+(--2)r). 

M.(x) := (34*‘ - 2^4+“) 

+ 2r(^T^ ~ ’ 

M3(x) := oDl-^cj){x) = I] Q {-^y{x - i)^+^ 

M4(x, /) := oDl-^Ml -x) = (^3(x - /)f' + - Of 

+ 2r(/3 + 3) (^^^^ ■ ^ - 9(x -; + 2)^+2 + 2(x -; + 3)^2) . 


kj^k—m ‘ 


a%‘ 


i!(-l)™ 
r(TO -I- (3) 


m+/3—1 


Then we have 


( 2 ‘^/ 2 (/)„^( 2 '^x)) = 2 '^( 3 / 2 -/ 3 )^^ , {i,a,) = ( 1 , a) or ( 2 , 6 ), 

oDl-^ ( 2 -^^^<t>{ 2 \ - fc)) = ^ 2 ‘^x - ft) , 

oDl-^ (2-^/^<j3b ( 2-^(1 - x))) = 2'^(3/2-/3)^^ ^ 2 '^a:, 2 ^^) . 

The similar formulae can also be derived for 2'^/^(/)a(2'^(x — 1)). 
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